Simulation

Load Grid

load('C:\Users\user\Documents\MATLAB\ProcejtGridVar.mat');
[~, ~, k] = gridLogicalIndices(G);
 
clf;
plotCellData(G, k, 'Edgealpha', 0.2);
plotFaces(G, find(G.faces.tag>0), 'Edgecolor', 'r', 'Facecolor', [0.8, 0.8, 0.8]);
colormap(0.7*jet + 0.3*ones(size(jet)));
view(-20, 50); axis off; camlight; material dull; camproj('perspective');

Define Petrophyisical Parameters

rng(1);
p = gaussianField(G.cartDims, [0.25, 0.35], 5, 2);
p = p(G.cells.indexMap);
 
K = logNormLayers(G.cartDims, [100, 400, 50]*milli*darcy, 'indices', [1, 3, 6, 10]);
K = K(G.cells.indexMap);
 
clf
 
subplot(2, 1, 1)
histogram(p);
xlabel('Porosity')
 
subplot(2, 1, 2)
histogram(convertTo(K, milli*darcy));
xlabel('Permeability [md]')

Porosity Distributions

clf;
% plotGrid(G, 'facecolor', 'none')
plotCellData(G, p);
title('Porosity distribution')
view(-20, 50);
colormap jet; colorbar('horiz')
axis off;

Permeability Distribution

% plotGrid(G, 'facecolor', 'none')
clf;
plotCellData(G, convertTo(K, milli*darcy));
title('Permeability distribution')
view(-20, 50)
colormap jet; colorbar('horiz');
axis off;
clf;
plotGrid(G, 'facecolor', 'none')
plotCellData(G, convertTo(K, milli*darcy), find(K>50*milli*darcy & p > 0.2));
title('Permeability distribution')
view(-20, 50)
colormap jet;
axis off;
clf;
plotGrid(G, 'facecolor', 'none')
plotCellData(G, convertTo(K, milli*darcy), find(K>50*milli*darcy & p > 0.2));
title('Permeability distribution')
 
colormap jet;
axis off;
 
view([-120.9 -26.0])
clf;
plotGrid(G, 'facecolor', 'none')
plotCellData(G, convertTo(K, milli*darcy), find(K>50*milli*darcy & p > 0.2));
title('Permeability distribution')
 
colormap jet;
axis off;
view([-216.4 5.8])

Define rock model

rock = makeRock(G, K, p);
 
cr = 1e-6/barsa;
p_r = 300*barsa;
pv_r = poreVolume(G, rock);
pv = @(p) pv_r .* exp( cr * (p - p_r) );
sv = @(p) G.cells.volumes - pv(p);

Define model for compressible fluid

mu0 = 5*centi*poise;
cmup = 2e-3/barsa;
cmut = 1e-3;
T_r = 300;
mu = @(p,T) mu0*(1+cmup*(p-p_r)).*exp(-cmut*(T-T_r));
 
beta = 1e-3/barsa;
alpha = 5e-3;
rho_r = 850*kilogram/meter^3;
rho_S = 750*kilogram/meter^3;
rho = @(p,T) rho_r .* (1+beta*(p-p_r)).*exp(-alpha*(T-T_r) );

Quantities for energy equation

Cp = 4e3;
Cr = 2*Cp;
Hf = @(p,T) Cp*(T)+(1-T_r*alpha).*(p-p_r)./rho(p,T);
Ef = @(p,T) Hf(p,T) - p./rho(p,T);
Er = @(T) Cr*T;
 

Wells

W = [];
 
% Constant Rate Wells
I = [1000, 900, 100, 450]*20/1000;
J = [900, 150, 900, 150]*20/1000;
R = [800, 600, 750, 430]*-meter^3/day;
nx = G.cartDims(1);
ny = G.cartDims(2);
nz = G.cartDims(3);
for i = 1: numel(R)
W = verticalWell(W, G, rock, I(i), J(i), 1:nz, 'name',['W', int2str(i)]);
end
 
NCFR = numel(R);
 
 
% Constant BHP Wells
 
clf;
plotCellData(G, k, 'Edgealpha', 0.2);
plotFaces(G, find(G.faces.tag>0), 'Edgecolor', 'r', 'Facecolor', [0.8, 0.8, 0.8]);
plotWell(G, W);
colormap(0.7*jet + 0.3*ones(size(jet)));
view(-20, 50); axis off; camlight; material dull; camproj('perspective');

Impose vertical equilibrium

gravity reset on; g = norm(gravity);
[z_0, z_max] = deal(0, max(G.cells.centroids(:,3)));
equil = ode23(@(z,p) g.* rho(p,T_r), [z_0, z_max], p_r);
p_init = reshape(deval(equil, G.cells.centroids(:,3)), [], 1); clear equil
T_init = ones(G.cells.num,1)*T_r;

Plot well and initial pressure

clf
plotCellData(G, p_init/barsa, 'EdgeColor', 'k');
title('Initial Pressure Distribution');
plotWell(G,W, 'height',10);
view(-20, 50), camproj perspective
colormap jet; colorbar('horiz');

Compute Transmissibilities

N = double(G.faces.neighbors);
intInx = all(N ~= 0, 2);
cf = G.cells.faces(:,1);
nf = G.faces.num;
N = N(intInx, :); % Interior neighbors
 
hT = computeTrans(G, rock); % Half-transmissibilities
Tp = 1 ./ accumarray(cf, 1 ./ hT, [nf, 1]); % Harmonic average
Tp = Tp(intInx); % Restricted to interior
 
kap = 4; % Heat conduction of granite
tmp = struct('perm',kap*ones(G.cells.num,1)); % Temporary rock object
hT = computeTrans(G, tmp); % Half-transmissibilities
Th = 1 ./ accumarray(cf, 1 ./ hT, [nf, 1]); % Harmonic average
Th = Th(intInx); % Restricted to interior
 

Define Discrete Operators

n = size(N,1);
C = sparse( [(1:n)'; (1:n)'], N, ones(n,1)*[-1 1], n, G.cells.num);
grad = @(x) C*x;
div = @(x) -C'*x;
avg = @(x) 0.5 * (x(N(:,1)) + x(N(:,2)));
upw = @(x,flag) x(N(:,1)).*double(flag)+x(N(:,2)).*double(~flag);

Define Flow Equations

% Writing in functional form means that v(p,T) is evaluated two
% times more than strictly needed if the whole definition of
% discrete equations written as a single function
gdz = grad(G.cells.centroids)*gravity()';
v = @(p,T) -(Tp./mu(avg(p),avg(T))).*(grad(p) - avg(rho(p,T)).*gdz);
pEq = @(p,T, p0,T0, dt) ...
(1/dt)*(pv(p).*rho(p,T) - pv(p0).*rho(p0,T0)) ...
+ div( avg(rho(p,T)).*v(p,T) );
hEq = @(p, T, p0, T0, dt) ...
(1/dt)*(pv(p ).*rho(p, T ).*Ef(p ,T ) + sv(p ).*Er(T ) ...
- pv(p0).*rho(p0,T0).*Ef(p0,T0) - sv(p0).*Er(T0)) ...
+ div( upw(Hf(p,T),v(p,T)>0).*avg(rho(p,T)).*v(p,T) ) ...
+ div( -Th.*grad(T));

Define Well Equations.

for i = 1:numel(W)
assert(numel(W(i).cells)==numel(unique(W(i).cells))); % each cell should only appear once
end
 
bhT = 200; % temperature of wells not used if production
 
p_conn = @(bhp,bhT, i) bhp + g*W(i).dZ.*rho(bhp,bhT); %connection pressures
q_conn = @(p,T, bhp, i) ...
W(i).WI .* (rho(p(W(i).cells),T(W(i).cells))./mu(p(W(i).cells),T(W(i).cells))) .* (p_conn(bhp,bhT, i)-p(W(i).cells));
 
 
rateEq = @(p,T, bhp, qS, i) qS-sum(q_conn(p, T, bhp, i))/rho_S;
BHPctrlEq = @(bhp, BHP) bhp-BHP;
RatectrlEq = @(qS, r) qS-r;
 
 

Initialize For Solution Loop

winx = zeros(numel(W), 1);
for j = 1:numel(W)
winx(j) = W(j).cells(1);
end
Qs_init = zeros(numel(W), 1);
 
[p_ad, T_ad, bhp_ad, qS_ad] = ...
initVariablesADI(p_init,T_init, p_init(winx), Qs_init);
nc = G.cells.num;
[pIx, TIx, bhpIx, qSIx] = deal(1:nc, nc+1:2*nc, ...
2*nc+1:2*nc+numel(W), 2*nc+numel(W)+1:2*nc+2*numel(W));
 
numSteps = 78; % number of time-steps
totTime = 365*day*1.5; % total simulation time
dt = totTime / numSteps; % constant time step
tol = 1e-5; % Newton tolerance
maxits = 10; % max number of Newton its
 
sol = repmat(struct('time',[], 'pressure',[], 'bhp',[], 'qS',[], ...
'T',[], 'qH',[]),[numSteps+1,1]);
sol(1) = struct('time', 0, 'pressure', value(p_ad), ...
'bhp', value(bhp_ad), 'qS', value(qS_ad), 'T', value(T_ad),'qH', 0);
 

Main loop

t = 0; step = 0;
hQ = zeros(numel(W), 1);
% hQ(i) = sum(value(hq));
hwb = waitbar(0,'Simulation..');
while t < totTime
t = t + dt;
step = step + 1;
fprintf('\nTime step %d: Time %.2f -> %.2f days\n', ...
step, convertTo(t - dt, day), convertTo(t, day));
 
% Newton loop
resNorm = 1e99;
p0 = value(p_ad); % Previous step pressure
T0 = value(T_ad);
nit = 0;
while (resNorm > tol) && (nit < maxits)
 
% Add source terms to homogeneous pressure equation:
eq1 = pEq(p_ad,T_ad, p0, T0,dt);
eq2 = hEq(p_ad,T_ad, p0, T0,dt);
for i = 1:numel(W)
qw = q_conn(p_ad,T_ad, bhp_ad(i), i);
eq1(W(i).cells) = eq1(W(i).cells) - qw;
hq = Hf(bhp_ad(i),bhT).*qw; %inflow not in this example
Hcells = Hf(p_ad,T_ad);
hq(qw<0) = Hcells(W(i).cells(qw<0)).*qw(qw<0);
hQ(i) = sum(value(hq));
eq2(W(i).cells) = eq2(W(i).cells) - hq;
Weq(i).rateEq = rateEq(p_ad,T_ad, bhp_ad(i), qS_ad(i), i);
if i <= NCFR
Weq(i).ctrlEq = RatectrlEq(qS_ad(i), R(i));
else
Weq(i).ctrlEq = BHPctrlEq(bhp_ad(i), bhps(i-NCFR));
end
end
% Collect all equations. Scale residual of energy equation
% Concatenate equations and solve for update:
 
eq = cat(eq1, eq2/Cp, Weq(:).rateEq, Weq(:).ctrlEq);
J = eq.jac{1}; % Jacobian
res = eq.val; % residual
upd = -(J \ res); % Newton update
 
% Update variables
p_ad.val = p_ad.val + upd(pIx);
T_ad.val = T_ad.val + upd(TIx);
bhp_ad.val = bhp_ad.val + upd(bhpIx);
qS_ad.val = qS_ad.val + upd(qSIx);
 
resNorm = norm(res);
nit = nit + 1;
fprintf(' Iteration %3d: Res = %.4e\n', nit, resNorm);
end
 
if nit > maxits
error('Newton solves did not converge')
else % store solution
sol(step+1) = struct('time', t, 'pressure', value(p_ad), ...
'bhp', value(bhp_ad), 'qS', value(qS_ad), ...
'T', value(T_ad), 'qH', hQ);
waitbar(t/totTime,hwb);
end
end
Time step 1: Time 0.00 -> 7.02 days
Iteration 1: Res = 3.6019e+02 Iteration 2: Res = 5.2653e+02 Iteration 3: Res = 1.9919e+02 Iteration 4: Res = 1.2039e+01 Iteration 5: Res = 4.2775e-02 Iteration 6: Res = 4.1497e-07
Time step 2: Time 7.02 -> 14.04 days
Iteration 1: Res = 2.2509e+02 Iteration 2: Res = 8.6686e-01 Iteration 3: Res = 1.4053e-05 Iteration 4: Res = 8.9263e-10
Time step 3: Time 14.04 -> 21.06 days
Iteration 1: Res = 1.6678e+02 Iteration 2: Res = 2.4684e-01 Iteration 3: Res = 1.0079e-06
Time step 4: Time 21.06 -> 28.08 days
Iteration 1: Res = 1.4810e+02 Iteration 2: Res = 1.3110e-01 Iteration 3: Res = 2.3178e-07
Time step 5: Time 28.08 -> 35.10 days
Iteration 1: Res = 1.3905e+02 Iteration 2: Res = 8.4086e-02 Iteration 3: Res = 7.6754e-08
Time step 6: Time 35.10 -> 42.12 days
Iteration 1: Res = 1.3432e+02 Iteration 2: Res = 6.0918e-02 Iteration 3: Res = 3.2365e-08
Time step 7: Time 42.12 -> 49.13 days
Iteration 1: Res = 1.3172e+02 Iteration 2: Res = 4.8609e-02 Iteration 3: Res = 1.7274e-08
Time step 8: Time 49.13 -> 56.15 days
Iteration 1: Res = 1.3022e+02 Iteration 2: Res = 4.1820e-02 Iteration 3: Res = 1.1765e-08
Time step 9: Time 56.15 -> 63.17 days
Iteration 1: Res = 1.2932e+02 Iteration 2: Res = 3.7963e-02 Iteration 3: Res = 9.6700e-09
Time step 10: Time 63.17 -> 70.19 days
Iteration 1: Res = 1.2874e+02 Iteration 2: Res = 3.5689e-02 Iteration 3: Res = 8.8410e-09
Time step 11: Time 70.19 -> 77.21 days
Iteration 1: Res = 1.2837e+02 Iteration 2: Res = 3.4285e-02 Iteration 3: Res = 8.3213e-09
Time step 12: Time 77.21 -> 84.23 days
Iteration 1: Res = 1.2811e+02 Iteration 2: Res = 3.3376e-02 Iteration 3: Res = 8.0968e-09
Time step 13: Time 84.23 -> 91.25 days
Iteration 1: Res = 1.2793e+02 Iteration 2: Res = 3.2760e-02 Iteration 3: Res = 7.9386e-09
Time step 14: Time 91.25 -> 98.27 days
Iteration 1: Res = 1.2781e+02 Iteration 2: Res = 3.2328e-02 Iteration 3: Res = 7.8217e-09
Time step 15: Time 98.27 -> 105.29 days
Iteration 1: Res = 1.2773e+02 Iteration 2: Res = 3.2019e-02 Iteration 3: Res = 7.7366e-09
Time step 16: Time 105.29 -> 112.31 days
Iteration 1: Res = 1.2767e+02 Iteration 2: Res = 3.1796e-02 Iteration 3: Res = 7.6913e-09
Time step 17: Time 112.31 -> 119.33 days
Iteration 1: Res = 1.2763e+02 Iteration 2: Res = 3.1636e-02 Iteration 3: Res = 7.6664e-09
Time step 18: Time 119.33 -> 126.35 days
Iteration 1: Res = 1.2761e+02 Iteration 2: Res = 3.1522e-02 Iteration 3: Res = 7.6764e-09
Time step 19: Time 126.35 -> 133.37 days
Iteration 1: Res = 1.2759e+02 Iteration 2: Res = 3.1446e-02 Iteration 3: Res = 7.6418e-09
Time step 20: Time 133.37 -> 140.38 days
Iteration 1: Res = 1.2759e+02 Iteration 2: Res = 3.1397e-02 Iteration 3: Res = 7.6524e-09
Time step 21: Time 140.38 -> 147.40 days
Iteration 1: Res = 1.2759e+02 Iteration 2: Res = 3.1372e-02 Iteration 3: Res = 7.6952e-09
Time step 22: Time 147.40 -> 154.42 days
Iteration 1: Res = 1.2759e+02 Iteration 2: Res = 3.1365e-02 Iteration 3: Res = 7.7291e-09
Time step 23: Time 154.42 -> 161.44 days
Iteration 1: Res = 1.2759e+02 Iteration 2: Res = 3.1372e-02 Iteration 3: Res = 7.7615e-09
Time step 24: Time 161.44 -> 168.46 days
Iteration 1: Res = 1.2759e+02 Iteration 2: Res = 3.1390e-02 Iteration 3: Res = 7.7940e-09
Time step 25: Time 168.46 -> 175.48 days
Iteration 1: Res = 1.2760e+02 Iteration 2: Res = 3.1418e-02 Iteration 3: Res = 1.1501e-08
Time step 26: Time 175.48 -> 182.50 days
Iteration 1: Res = 1.2760e+02 Iteration 2: Res = 3.1453e-02 Iteration 3: Res = 7.8678e-09
Time step 27: Time 182.50 -> 189.52 days
Iteration 1: Res = 1.2760e+02 Iteration 2: Res = 3.1494e-02 Iteration 3: Res = 7.9286e-09
Time step 28: Time 189.52 -> 196.54 days
Iteration 1: Res = 1.2761e+02 Iteration 2: Res = 3.1540e-02 Iteration 3: Res = 7.9607e-09
Time step 29: Time 196.54 -> 203.56 days
Iteration 1: Res = 1.2761e+02 Iteration 2: Res = 3.1589e-02 Iteration 3: Res = 8.0005e-09
Time step 30: Time 203.56 -> 210.58 days
Iteration 1: Res = 1.2761e+02 Iteration 2: Res = 3.1641e-02 Iteration 3: Res = 8.0729e-09
Time step 31: Time 210.58 -> 217.60 days
Iteration 1: Res = 1.2760e+02 Iteration 2: Res = 3.1695e-02 Iteration 3: Res = 8.1110e-09
Time step 32: Time 217.60 -> 224.62 days
Iteration 1: Res = 1.2760e+02 Iteration 2: Res = 3.1751e-02 Iteration 3: Res = 8.1990e-09
Time step 33: Time 224.62 -> 231.63 days
Iteration 1: Res = 1.2759e+02 Iteration 2: Res = 3.1808e-02 Iteration 3: Res = 8.2558e-09
Time step 34: Time 231.63 -> 238.65 days
Iteration 1: Res = 1.2759e+02 Iteration 2: Res = 3.1866e-02 Iteration 3: Res = 8.2849e-09
Time step 35: Time 238.65 -> 245.67 days
Iteration 1: Res = 1.2758e+02 Iteration 2: Res = 3.1924e-02 Iteration 3: Res = 8.3921e-09
Time step 36: Time 245.67 -> 252.69 days
Iteration 1: Res = 1.2757e+02 Iteration 2: Res = 3.1983e-02 Iteration 3: Res = 8.4257e-09
Time step 37: Time 252.69 -> 259.71 days
Iteration 1: Res = 1.2755e+02 Iteration 2: Res = 3.2043e-02 Iteration 3: Res = 8.4975e-09
Time step 38: Time 259.71 -> 266.73 days
Iteration 1: Res = 1.2754e+02 Iteration 2: Res = 3.2103e-02 Iteration 3: Res = 8.5546e-09
Time step 39: Time 266.73 -> 273.75 days
Iteration 1: Res = 1.2753e+02 Iteration 2: Res = 3.2163e-02 Iteration 3: Res = 8.5914e-09
Time step 40: Time 273.75 -> 280.77 days
Iteration 1: Res = 1.2751e+02 Iteration 2: Res = 3.2223e-02 Iteration 3: Res = 8.6912e-09
Time step 41: Time 280.77 -> 287.79 days
Iteration 1: Res = 1.2750e+02 Iteration 2: Res = 3.2283e-02 Iteration 3: Res = 8.7229e-09
Time step 42: Time 287.79 -> 294.81 days
Iteration 1: Res = 1.2748e+02 Iteration 2: Res = 3.2343e-02 Iteration 3: Res = 8.8178e-09
Time step 43: Time 294.81 -> 301.83 days
Iteration 1: Res = 1.2746e+02 Iteration 2: Res = 3.2403e-02 Iteration 3: Res = 8.8617e-09
Time step 44: Time 301.83 -> 308.85 days
Iteration 1: Res = 1.2744e+02 Iteration 2: Res = 3.2464e-02 Iteration 3: Res = 8.9523e-09
Time step 45: Time 308.85 -> 315.87 days
Iteration 1: Res = 1.2742e+02 Iteration 2: Res = 3.2524e-02 Iteration 3: Res = 9.0143e-09
Time step 46: Time 315.87 -> 322.88 days
Iteration 1: Res = 1.2740e+02 Iteration 2: Res = 3.2585e-02 Iteration 3: Res = 9.0862e-09
Time step 47: Time 322.88 -> 329.90 days
Iteration 1: Res = 1.2738e+02 Iteration 2: Res = 3.2645e-02 Iteration 3: Res = 9.1276e-09
Time step 48: Time 329.90 -> 336.92 days
Iteration 1: Res = 1.2736e+02 Iteration 2: Res = 3.2706e-02 Iteration 3: Res = 9.1923e-09
Time step 49: Time 336.92 -> 343.94 days
Iteration 1: Res = 1.2734e+02 Iteration 2: Res = 3.2766e-02 Iteration 3: Res = 9.2713e-09
Time step 50: Time 343.94 -> 350.96 days
Iteration 1: Res = 1.2731e+02 Iteration 2: Res = 3.2827e-02 Iteration 3: Res = 9.3446e-09
Time step 51: Time 350.96 -> 357.98 days
Iteration 1: Res = 1.2729e+02 Iteration 2: Res = 3.2888e-02 Iteration 3: Res = 9.3990e-09
Time step 52: Time 357.98 -> 365.00 days
Iteration 1: Res = 1.2727e+02 Iteration 2: Res = 3.2949e-02 Iteration 3: Res = 9.4887e-09
Time step 53: Time 365.00 -> 372.02 days
Iteration 1: Res = 1.2724e+02 Iteration 2: Res = 3.3010e-02 Iteration 3: Res = 9.5509e-09
Time step 54: Time 372.02 -> 379.04 days
Iteration 1: Res = 1.2722e+02 Iteration 2: Res = 3.3072e-02 Iteration 3: Res = 9.6256e-09
Time step 55: Time 379.04 -> 386.06 days
Iteration 1: Res = 1.2719e+02 Iteration 2: Res = 3.3134e-02 Iteration 3: Res = 9.6960e-09
Time step 56: Time 386.06 -> 393.08 days
Iteration 1: Res = 1.2717e+02 Iteration 2: Res = 3.3196e-02 Iteration 3: Res = 9.7706e-09
Time step 57: Time 393.08 -> 400.10 days
Iteration 1: Res = 1.2714e+02 Iteration 2: Res = 3.3258e-02 Iteration 3: Res = 9.8651e-09
Time step 58: Time 400.10 -> 407.12 days
Iteration 1: Res = 1.2712e+02 Iteration 2: Res = 3.3320e-02 Iteration 3: Res = 9.9306e-09
Time step 59: Time 407.12 -> 414.13 days
Iteration 1: Res = 1.2709e+02 Iteration 2: Res = 3.3383e-02 Iteration 3: Res = 9.9884e-09
Time step 60: Time 414.13 -> 421.15 days
Iteration 1: Res = 1.2706e+02 Iteration 2: Res = 3.3446e-02 Iteration 3: Res = 1.0065e-08
Time step 61: Time 421.15 -> 428.17 days
Iteration 1: Res = 1.2704e+02 Iteration 2: Res = 3.3510e-02 Iteration 3: Res = 1.0161e-08
Time step 62: Time 428.17 -> 435.19 days
Iteration 1: Res = 1.2701e+02 Iteration 2: Res = 3.3574e-02 Iteration 3: Res = 1.0235e-08
Time step 63: Time 435.19 -> 442.21 days
Iteration 1: Res = 1.2698e+02 Iteration 2: Res = 3.3638e-02 Iteration 3: Res = 1.0300e-08
Time step 64: Time 442.21 -> 449.23 days
Iteration 1: Res = 1.2696e+02 Iteration 2: Res = 3.3702e-02 Iteration 3: Res = 1.0382e-08
Time step 65: Time 449.23 -> 456.25 days
Iteration 1: Res = 1.2693e+02 Iteration 2: Res = 3.3767e-02 Iteration 3: Res = 1.0466e-08
Time step 66: Time 456.25 -> 463.27 days
Iteration 1: Res = 1.2690e+02 Iteration 2: Res = 3.3833e-02 Iteration 3: Res = 1.0556e-08
Time step 67: Time 463.27 -> 470.29 days
Iteration 1: Res = 1.2688e+02 Iteration 2: Res = 3.3899e-02 Iteration 3: Res = 1.0643e-08
Time step 68: Time 470.29 -> 477.31 days
Iteration 1: Res = 1.2685e+02 Iteration 2: Res = 3.3965e-02 Iteration 3: Res = 1.0711e-08
Time step 69: Time 477.31 -> 484.33 days
Iteration 1: Res = 1.2682e+02 Iteration 2: Res = 3.4032e-02 Iteration 3: Res = 1.0780e-08
Time step 70: Time 484.33 -> 491.35 days
Iteration 1: Res = 1.2679e+02 Iteration 2: Res = 3.4099e-02 Iteration 3: Res = 1.0873e-08
Time step 71: Time 491.35 -> 498.37 days
Iteration 1: Res = 1.2676e+02 Iteration 2: Res = 3.4167e-02 Iteration 3: Res = 1.0962e-08
Time step 72: Time 498.37 -> 505.38 days
Iteration 1: Res = 1.2674e+02 Iteration 2: Res = 3.4235e-02 Iteration 3: Res = 1.1059e-08
Time step 73: Time 505.38 -> 512.40 days
Iteration 1: Res = 1.2671e+02 Iteration 2: Res = 3.4304e-02 Iteration 3: Res = 1.1150e-08
Time step 74: Time 512.40 -> 519.42 days
Iteration 1: Res = 1.2668e+02 Iteration 2: Res = 3.4373e-02 Iteration 3: Res = 1.1235e-08
Time step 75: Time 519.42 -> 526.44 days
Iteration 1: Res = 1.2665e+02 Iteration 2: Res = 3.4443e-02 Iteration 3: Res = 1.1328e-08
Time step 76: Time 526.44 -> 533.46 days
Iteration 1: Res = 1.2662e+02 Iteration 2: Res = 3.4513e-02 Iteration 3: Res = 1.1425e-08
Time step 77: Time 533.46 -> 540.48 days
Iteration 1: Res = 1.2660e+02 Iteration 2: Res = 3.4585e-02 Iteration 3: Res = 1.1502e-08
Time step 78: Time 540.48 -> 547.50 days
Iteration 1: Res = 1.2657e+02 Iteration 2: Res = 3.4656e-02 Iteration 3: Res = 1.1605e-08
Time step 79: Time 547.50 -> 554.52 days
Iteration 1: Res = 1.2654e+02 Iteration 2: Res = 3.4728e-02 Iteration 3: Res = 1.1698e-08
close(hwb)

Visualization

Pressure Evolution

clf;
steps = [2 5 10 25];
%steps = floor(1:(numel(sol)/6-eps):numel(sol));
p = vertcat(sol(:).pressure);
cax=[min(p) max(p)]./barsa;
for i=1:4
subplot(2,2,i);
set(gca,'Clipping','off');
plotCellData(G, sol(steps(i)).pressure/barsa, 'EdgeColor',.5*[1 1 1]);
plotWell(G, W, 'FontSize',12);
view(70,45), camproj perspective
caxis(cax);
axis tight off; zoom(1.1)
text(400,170,-8,[num2str(round(steps(i)*dt/day)) ' days'],'FontSize',12);
end
sgtitle('Pressure Evolution')
colorbar('horiz','Position',[.1 .05 .8 .025]);
colormap(jet(55));
clf
steps = [30 45 60 numel(sol)];
%steps = floor(1:(numel(sol)/6-eps):numel(sol));
p = vertcat(sol(:).pressure);
cax=[min(p) max(p)]./barsa;
for i=1:4
subplot(2,2,i);
set(gca,'Clipping','off');
plotCellData(G, sol(steps(i)).pressure/barsa, 'EdgeColor',.5*[1 1 1]);
plotWell(G, W, 'FontSize',12);
view(70,45), camproj perspective
caxis(cax);
axis tight off; zoom(1.1)
text(400,170,-8,[num2str(round(steps(i)*dt/day)) ' days'],'FontSize',12);
end
sgtitle('Pressure Evolution')
colorbar('horiz','Position',[.1 .05 .8 .025]);
colormap(jet(55));
Res.gif

Temperature Evolution

clf;
steps = [2 5 10 numel(sol)];
T = vertcat(sol(:).T); %-T_r;
cax=[min(T) max(T)];
for i=1:numel(steps)
subplot(2,2,i);
set(gca,'Clipping','off');
plotCellData(G, sol(steps(i)).T, 'EdgeColor',.5*[1 1 1]);
plotWell(G, W, 'FontSize',12);
view(70,45), camproj perspective
caxis(cax);
axis tight off; zoom(1.1)
text(400,170,-8,[num2str(round(steps(i)*dt/day)) ' days'],'FontSize',12);
end
sgtitle('Temperature Evolution')
h=colorbar('horiz','Position',[.1 .05 .8 .025]);
TemperatureEvolution.gif

Plot production rate and bhp

clf;
set(gca,'FontSize',20);
 
rate = abs([sol(2:end).qS]*day);
BHPs = [sol(2:end).bhp]/barsa;
for i = 1:numel(W)
subplot(2, 2, i);
stairs([sol(2:end).time]/day, rate(i, :),'LineWidth',2);
ylabel('Rate [meter3/day]')
yyaxis right;
stairs([sol(2:end).time]/day, BHPs(i, :),'LineWidth',2);
title(W(i).name);
xlabel('Time [days]');
ylabel('BHP [bars]')
grid on;
 
end

Plot Cumulative Production

clf;
cumprod = zeros(1 , numel(sol) - 1);
for i = 1:numel(W)
cumprod = cumprod + cumtrapz([sol(2:end).time]/day, rate(i, :));
end
stairs([sol(2:end).time]/day,cumprod,'LineWidth',2);
title('Cumpulative Production');
xlabel('Time [days]')
ylabel('meter3')
grid on;

Plot Pressures and Temperatures

i = 1;2;
ns = numel(sol);
Tw = nan(ns,numel(W(i).cells));
Pw = nan(ns,numel(W(i).cells));
[t,Pm,PM,Pa,Tm,TM,Ta] = deal(nan(ns,1));
i = 1;
for j=1:numel(sol)
t(j) = sol(j).time/day;
Tw(j,:)=sol(j).T(W(i).cells);
Tm(j) = min(sol(j).T);
TM(j) = max(sol(j).T);
Ta(j) = mean(sol(j).T);
Pw(j,:)=sol(j).pressure(W(i).cells)./barsa;
Pm(j) = min(sol(j).pressure./barsa);
PM(j) = max(sol(j).pressure./barsa);
Pa(j) = mean(sol(j).pressure./barsa);
end
clf;
plot(t,Pm,t,Pa,t,PM,t,Pw,'.k','LineWidth',2);
legend('min(p)', 'avg(p)', 'max(p)', 'wells');
title('Pressure');
xlabel('Time [days]')
ylabel('bars')
 
 
grid on;

Compute the three different expansion temperatures

[p,T] = initVariablesADI(p_r,T_r);
dp = Pm(end)*barsa-p_r;
 
% Joule-Thomson
hf = Hf(p,T);
dHdp = hf.jac{1};
dHdT = hf.jac{2};
Tjt = T_r - dHdp*dp/dHdT;
hold on, plot(t([1 end]), [Tjt Tjt],'--k'); hold off
text(t(5), Tjt+.25, 'Joule-Tompson');
 
% linearized adiabatic temperature
hf = Ef(p,T) + value(p)./rho(p,T);
dHdp = hf.jac{1};
dHdT = hf.jac{2};
Tab = T_r - dHdp*dp/dHdT;
hold on, plot(t([1 end]), [Tab Tab],'--k'); hold off
text(t(2), Tab-.35, 'Adiabatic expansion');
 
% free expansion
hf = Ef(p,T);
dHdp = hf.jac{1};
dHdT = hf.jac{2};
Tfr = T_r - dHdp*dp/dHdT;
hold on, plot(t([1 end]), [Tfr Tfr],'--k'); hold off
text(t(end/2), Tfr+.25, 'Free expansion');
set(gca,'XLim',t([1 end]));